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I. INTRODUCTION 



Gravitational wave astronomy is an active field of research which promises to open up a new window into 
the physical universe Q], 0], The current and future laser interferometric gravitational 

wave detectors, high precision pulsar timing, along with measurements of the anisotropics in the temperature 
and polarization of the Cosmic Microwave Background have the potential to discover gravitational waves in 
a broad range of frequencies in the near future (see [3| for recent discussions). 

In this paper we shall be mainly interested in pulsar timing as a laboratory for gravitational wave physics. 
Propagation of pulsar signal through space-time perturbed by gravitational waves results in appearance of 
anomalous timing residuals (i.e. differences between observed and theoretically predicted times of arrival). 
Pulsar timing provides a unique tool for observing gravitational waves in low-frequency band (10~^ Hz < 
fgw < 10~^ Hz) [lli, [l^, [l^, [l^, [l^, [l^, The main sources of gravitational waves at these frequencies 
are expected to be of extragalactic orig in. The strongest sources would be supermassivc black hole binaries 
in the center of galaxies [l3], [la], 0|- Relic gravitational waves, which are the remnants from the 
early history of the universe, may also contribute a significant fraction to the gravitational wave background 



at these frequencies 



2l| , [22l | . Pulsar timing could also measure gravitational waves from superstrings [23| , 



as well as several other exotic sources 



24|. 



The main methods to detect gravitational waves are based on the analysis of their interaction with 



electromagnetic fields 



25 



26[, 27|. The interaction of gravitational waves with electromagnetic waves 



leaves measurable imprints on the latter. For example, the phase variations in the electromagnetic wave 
propagatin g in the field of a gravitational wave, and its implications for space radio intcrferometry were 



studied in 



28l | (see also analyzing these phase variations in a situation when the speed of 



gravitational waves could be smaller than the speed of light, the authors introduced the concept of "surfing 
effect" and studied its implications for the precision intcrferometry measurements. In this paper wc shall 
consider the implications of the surfing effect for pulsar timing measurements. As we shall show, due to the 
transverse nature of gravitational waves, the surfing effect can lead to enormous observable pulsar timing 
residuals if the speed of gravitational waves is smaller than the speed of electromagnetic waves. Wc shall use 
this fact, along with the expected precision of pulsar timing measurements, to place stringent upper limits 
on the parameter e = (c — Vgw ) /c which characterizes the deviation of speed of gravitational waves from the 
speed of light. We show that, for a realistic gravitational wave background and a reasonable time duration 
of observations, the achievable limits are e < 0.4%. Constraining the speed of gravitational waves is an 
interesting experimental challenge attracting much theoretical and experimental interest 

[3, H. H- We 

argue that the constraint on e from pulsar timing would provide the strongest current limitations on the 
deviation of speed of gravitational waves from speed of light. 

It is worth mentioning that the surfing effect considered in this paper is quite generic. The surfing 
effect occurs in any physical situation where the phase speed of gravitational waves is smaller than the 
phase speed of electromagnetic waves jsol l. For example, this is the case in theories which predict 
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a non vanishing rest mass for graviton 



34| . 3ll |. 35| . Although, generically, these theories predict extra 



polarization states for gravitational waves, in our work we shall restrict our analysis to effects caused only 
by transverse tracelcss (TT) gravitational waves. Another possible scenario for the surfing effect to arise is 
to consider the interaction of gravitational waves and electromagnetic waves in the presence of plasma. In 

Dccd of light 



this case the phase speed of gravitational waves remains unchanged and is equal to c (i.e. the sr 
in vacuum), while the phase speed of electromagnetic waves becomes generally greater than c 

The plan of the paper is as follows. We shall begin in Section [IT] with the analysis of propagation of an 
electromagnetic wave in the field of a single monochromatic plane gravitational wave. Wc shall calculate the 
timing residuals due to a single gravitational wave and discuss the manifestations and physical consequences 
of the surfing effect. In Section [llll we generalize the surfing effect for the case of an arbitrary gravitational 
wave field. We derive the statistical properties of the timing residual signal based on the statistical properties 
of the gravitational wave field. In Section HVl we calculate the achievable constraints on e depending on the 
strength of the gravitational wave background characterized by energy density parameter In Section 

|V|we study the physical consequences of the surfing effect in pulsar timing. We show that the gravitational 
wave background from extragalactic black holes allows to place strong limits on e. Furthermore we show 
that the surfing effect can also place a strong upper bound on the mass of graviton. Finally, we conclude 
the paper in Section IVII with a summary of the main results of this work. 



II. PULSAR TIMING RESIDUALS FOR A SINGLE MONOCHROMATIC GRAVITA- 
TIONAL WAVE 



In this paper wc shall be working in the framework of a slightly perturbed Minkowski space time with 
coordinates = {ct, .t') and the metric given by 



-c^dt^ + (% + hij) dx^dx^ , 



(1) 



where hij is the gravitational wave perturbation. For clarity and in order to gain physical insight into the 
problem, in this section, wc shall consider the case of a single monochromatic plane gravitational wave. In 
the next section, we shall generalize our analysis to the case of an arbitrary gravitational wave field. For a 
monochromatic gravitational wave the metric perturbation hij takes the form 



25|, [26| 



h pijC 



(2) 



where h is the amplitude of the gravitational wave, k^^ ~ (fco, fci) is the wave vector, and pik is the polarization 
tensor of the gravitational wave. Introducing a set of two mutually orthogonal unit vectors U and rrii 
orthogonal to the wave vector ki , the polarization tensor pik has the form 25| , 26 1 



Pik = -[k± rui) {Ik ± rrik) , 



(3) 
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where ± corresponds to the two independent states of chxular polarization. Due to the transverse and 
traceless nature of gravitational waves, the polarization tensor satisfies the following conditions 

P^kk' = 0, m<5''- = 0. (4) 

1 /2 

For further discussion, it is convenient to introduce the wavenumber k — (^Sijk^y^ , and a unit vector in 
the direction of wave propagation k' ~ k'' /k. The wavelength of the gravitational wave is related to the 
wavenumber by the equality k = 2TT/Xgw. The frequency of the gravitational wave fgw is related to the time 
component of the wave vector through the relation ko = 2TTfgm/c. 

The speed of a gravitational wave is determined by relationship Vg^ = fgw^gw In General Relativity 
gravitational waves travel at the speed of light, i.e. Vg^ = c, which implies a relationship (dispersion rela- 
tionship) k — ko- In order to analyze the possibility Vg^ ^ c, let us introduce a phcnomcnological parameter 
e describing the relative deviation of Vgw from speed of light c 

e = ^ — where Vgn, = fgy^Xgy, = ^ = c (1 - e) . (5) 

The quantity e has been introduced as a phenomenological parameter, and thus the analysis that follows is 

valid for any theory that predicts gravitational waves with Vg^ ^ c. Particularly, of interest are modifications 

of General Relativity that predict massive gravitons. For these models, e can be related to the rest mass of 

the graviton rUg through the relation 

^_ fickp ^ TUgC _ nigC^ 

hcko + nigC^ hkf) 27Thfgiu 

Let us move our attention to pulsar timing measurements. The effect of a gravitational wave upon the 



measured frequency of pulsar signal is given by , 



D 

(7) 



vq 2c J \ dt 



path 



where vq is the unperturbed pulsar frequency in the absence of gravitational waves and Ai'{t) ~ ^{t) — vq 
is the variation of pulsar frequency due to the presence of a gravitational wave. D is the distance from the 
pulsar to the observer, integration variable s is the distance parameter along the unperturbed light ray path 
from pulsar to the observer, is the unit vector tangent along this path (i.e. unit vector in the direction 
from pulsar to the observer), and the subscript indicates the integration along this path. The unperturbed 
light ray path is given by 

t(s) = i--, x'{s) =x' -e's, (8) 
c 

where t and x'' determine the time and position of the observation. Without loss of generality we can set 
= by choosing a spatial coordinate system with observer at its origin. 

Substituting the path ([8]) into (O, taking into account ([2|) and (O, after straight forward integration we 
arrive at 

_ gi(l-e-fcie')fcD 



vq 2 



1 - e - heJ 



(9) 
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The pulsar timing measurements customarily measure the timing residuals, i.e. the difference between 
the actual pulse arrival times and times predicted from a spin-down model for a pulsar [l^ . [l^. The 
variations in the measured frequency, due to the presence of a gravitational wave, will cause an anomalous 
timing residual R{t) in the pulse arrival time given by 12 1 

t 



R{t) 



dt 



t-T 



^0 



(10) 



where T is the time of observations, and the residual R{t) is measured in seconds. Substituting expression 
into (jlOp . we get for the timing residual due to a single monochromatic gravitational wave 



R[t) = e'e^py e-**^(i-^)^* f 1 - e*Mi-OcT\ 

2kc V / 



(11) 



Before proceeding further, let us analyze the above expression. The expression in the square brackets on 
the right side of (fTTI) becomes large (proportional to kD ^ D/Xg^) when ( 1 — e — fc.; 



0, i. e. 

R{t) w e^e^py g-»fc(i-E)ct _ g»fc(i-^)cT^ (l + O (,5))] , for ^ = (^1 - e - fc,e') kD < 1. (12) 

Hence, for gravitational waves traveling in a direction at a sufficiently small angle to the direction from the 
pulsar, i.e. fciC* w (1 — e), there is a resonance inreasc in the expression for timing residual. In the case 
when e = this does not lead to a growth of the timing residual R{t) itself, due to the transverse nature 
of the gravitational wave (since e''eJpij — > when fc^e* — > 1. see expression (HU). On the other hand, if 
e 7^ 0, the expression for R{t) increases significantly for kic'' « (1 — e). The resonance occurs when the 
signal from the pulsar "surfs" along the gravitational wave, i.e. travels at a small angle cos 6* » (1 — e) to 



30| we call 



the gravitational wave. This picture is reminiscent of wave surfing, so for this reason following 
this effect, of a resonant increase va.R{t), as the surfing effect. It is worth noticing that the above analysis 



closely reseixibles considerations in 



30|, where the surfing effect manifested itself in the resonance growth 



of the phase variation of electromagnetic waves, leading to an observable angular displacement of distant 
quasars. In the current work, we are analyzing the signature of the surfing effect in pulsar timing residuals. 



III. PULSAR TIMING RESIDUALS FOR AN ARBITRARY GRAVITATIONAL WAVE 
FIELD 

In the previous section we calculated the timing residual due to a single plane monochromatic gravi- 
tational wave. In this section we shall generalize our analysis to an arbitrary gravitational wave field. In 
general, an arbitrary gravitational wave field can be decomposed into spatial Fourier modes 



hij 



{t,x') - / d^k [hs{k\t) (fc')e'''='"' + K{k\t) pI (fc')e-'^'"'] , (13) 

s=l,2 
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where o?^k denotes the integration over all possible wave vectors, and s = 1, 2 corresponds to the two linearly 
independent modes of polarization satisfying the orthogonality condition 



PijP =Ss 



(14) 



The mode function hs{k^,t) correspond to plane monochromatic waves 

hs{k\t) = hs{k') e^''''-^-'^"'' (15) 

Due to the linear nature of the problem, following the decomposition (|13[) . the total timing residual due to 
an arbitrary gravitational wave field, can be presented in the following manner 



R{t)= / rf^k ^ \^h,{k')R{t;k\s) + hl{k')R*{t]k\s) 



s=l,2 



(16) 



Using the results of the previous section, the contribution from a single Fourier component R{t;k'^,s) is 
given by 



2kc 



l-e- k, 



(17) 



where the tilde over R in the above expression is introduced to indicate explicit factoring out of the gravi- 
tational wave amplitude h compared with 

In general, if we have the information about the mode functions hs{k^), using expressions (|16p and (|17p 
we can calculate the expected timing residual for an arbitrary gravitational wave field. In most of the 
practically interesting cases we do not have such a complete knowledge of the gravitational wave field, but 
are restricted to the knowledge of its statistical properties. To proceed, let us assume the following statistical 
properties 



< h,{k') >= 0, < h,{k') hl,{k'') iM^5,,,^3(fc'' _ fc'O, 



(18) 



167rfc-^ 

where the brackets denote ensemble averaging over all possible realizations, and Ph{k) is the metric power 
spectrum per logarithmic interval of k. These conditions correspond to a stationary statistically homoge- 
neous and isotropic gravitational wave field. 

The positing of the statistical properties of the gravitational wave field p8|l allows us to calculate the 
statistical properties of the timing residual R{t). Using (|16p and p^ . and taking into account the orthogo- 
nality property (|14p . after straight forward calculations, we arrive at the following statistical properties for 
the timing residual R 



< R{t) > = 0, 



dk 



<R^{t) > = I —P,,ik)R^{k), 



where we have introduced the transfer function 



(19a) 
(19b) 

(20) 
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In the above expression dfi represents integration over the possible directions of gravitational wave (i.e. d'^k = 
k^dkdfl). From (fT7|) and ([^0]) it follows that the transfer function E?{k) does not depend on time variable 
t, which is a reflection of the stationarity of the underlying gravitational wave field. 

The expression for the transfer function can be explicitly calculated. In order to do this, let us firstly 
introduce a spherical coordinate system (6', (/>) related to the spatial coordinates {x'} (following notations 
of 37[). Without loss of generality, we can assume that our spatial coordinate system is chosen such 



that the unit vector from the pulsar to the observer points in the north-pole direction, i.e. e* = (0,0,1). 
Let us also introduce the quantity /i = cos 9 = Eik^ , characterizing the angle between the direction of 
gravitational wave propagation and the direction from pulsar to the observer. Furthermore, let cf) denote 
the azimuthal angle that is subtended by A:* projected onto the (a;^, a:;^)-plane, i.e. k^ ~ cos(/)sin0 and 
k/^ — sin sin 0. Introducing and ef which arc the meridian and azimuthal unit vectors perpendicular to 
the gravitational wave wavevector ki respectively, the polarization tensors for gravitational waves ([3]) take 
the form (fc') = (e^ ± ief )(e^ ± ze^)/2, with ± corresponding to the two i ndep endent circularly polarized 
degrees of freedom s = 1,2 (for a detailed discussion see for example 38[, 39[). Taking into account the 
relation 

eVP.,= i(l-/.^)e±^^^ (21) 

substituting (jl7p into (j20p and setting dVL ~ dfid(j), after integration over 0, we arrive at the expression for 
the transfer function 

(22) 



1 2 f kcT 
2fcV"" [ — 



The integrand under the integral in the above expression is illustrated in Figure [TJ As can been seen, when 
e 7^ 0, the predominant contribution to the integral comes from the resonance region /i w (1 — e). Thus, in 
this case, the predominant contribution to the timing residual < > comes from "surfing" gravitational 
waves, i.e. waves for which /i « (1 — e). In the physically interesting limit e — > and kD —^oowe can 
calculate the integral in (22) explicitly. We refer the reader to Appendix [XI for details of this calculation. 
The result is as follows 





3 , 




1 + -TTE^kD 





(23) 



The above expression allows us to simply quantify the condition for the surfing effect to be dominant, 
e^kD 3> 1. As we shall show in the next section, given the precision level of the current and planned pulsar 
timing measurements, the surfing effect allows to place significant constraints on the e parameter. 

Before proceeding, it is instructive to compare the results of this section with the results of [s^ . More 
specifically, it is interesting to compare expression (|23p for the transfer function of timing residuals with its 
counterpart expression (29) in [s^ for the transfer function of angular displacement 



^a^{k) w — [1 + 5Tre^ kD] 
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\J.I u.o 

|j. = cose = cose 



FIG. 1: The illustration of the resonance effect, present for e 7^ 0. The graphs show integrand in expression 
(|22p . For the case e 7^ the integrand sharply peaks at angle /i w {1 — e) (solid red line), while for the case 
e = the effect is absent (dashed blue line). In the case of e 7^ 0, the gravitational waves travelling at an 
angle cos 6' « (1 — e) to the line of sight are the predominant contributors to the surfing effect. The figure 
on the left shows the integrand for the whole region of /i, while the figure on the right zooms into the region 
around the resonance. 



Apart from the differing factors in front of the square brackets in the two expression, the crucial difference 
is the differing powers of e. In the present work, the surfing effect manifests in the term e^kD in the 
square brackets of In [30|, the surfing effect manifests in the term e^kD term in the square brackets 

of (29). The extra factor of e arose due to the geometrical specificity of interferometric observations of 
phase difference at the ends the interferometric system (see [s^ for details). The main consequences of this 
difference are twofold. Firstly, equivalent constraints on e require smaller distance to the source in the case 
of pulsar timing compared with interferometric observations. This is reflected in the fact that in the present 
work we focus on galactic pulsars, where as js^ focused on high redshift quasars. Secondly, the condition 
for surfing effect to dominate is different in the two contexts. This condition, characterized by the value of 
e, (see expression below and expression (32) in [s^), places the lower limit on the potentially possible 
bounds on e. This limiting bound is lower for interferometry measurements (e» ~ 2.3 x 10^'*) than for pulsar 
timing measurements (e» ~ 3.2 x 10"^). Even so. due to exceptional precision, the experimentally achievable 
bounds on e from pulsar timing measurements would be more stringent. 



IV. UPPER LIMITS ON THE SPEED OF GRAVITATIONAL WAVES 

Let us now turn our attention to the various cosmological and astrophysical candidates for a stochastic 
gravitational wave background and their contribution to the surfing effect in pulsar timing measurement. 
Analyzing their magnitude, we shall study the achievable upper limits on e that these backgrounds could 
place. 
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The stochastic gravitational wave field may be characterized by the dimcnsionless strain amplitude hc{f) 
which is related to the power spectrum in the following way 

ck 



hciD^^Kik), where / = —(!-£). (24) 

The quantity hd f) is the root-mean value of the gravitational wave amplitude in a unit logarithmic interval 
of frequencies. For analyzing the stochastic gravitational wave fields, it is also customary to introduce the 
density parameter Qg^ to characterize the strength of the gravitational wave field Q , Q , 3 • ^gw is related 
to the power spectrum Ph{k) and strain hc{f) by the relation 



= = m) (25) 

where kn = 27r//f/c = 2ttHo/c, and Ho is the current Hubble parameter. The density parameter Qg^ 
is the current day ratio of energy density of gravitational waves (per unit logarithmic interval in k) to 
the critical density of the Universe Pcrit = Sc^Ho/SttG. Below, for numerical estimations, we set Hubble 
parameter Ho = 75 ^ /Mpc. Note, that the above definition is valid for stationary gravitational wave 
backgrounds. In cosmological context, when considering relic gravitational, waves this definition modifies 
to flgwik) — ^ [j^l Ph{k) due to the non-stationary (standing wave) nature of relic gravitational waves 



3 \r^TFf / 

(see for example [3[). 

For simplicity, in the numerical estimations below, we shall assume a simple power law behaviour for he 
which is equivalent to a power law spectrum for the density parameter fig^, 

hc{f)^ho{fo){j-^ , f^3.(fc)-17<,„(fc„)(^) (26) 



where 

^^9»(fco) = ^ (^^j hiifo), ko^^, nT = 2(1 + a). (27) 

Although restricted, this form of spectrum is a good approximation for a large variety of models in gravita- 
tional wave frequency range of our interest. For example, this type of power law spectrum, with a = —2/3, is 
produced by the extragalactic coalescing super massive binary black hole systems \v\ . In cosmological con- 
text, this type of a power spectrum, with spectral index a at the current epoch, arises due to the evolution of 
relic gravitational waves with a primordial spectral index equal to 2(1 + a), (i.e. Ph{k)\prim oc fc2^^+")) 2l| . 
The flat, scale invariant power spectrum (also known as Harrison- Zeldovich power spectrum) corresponds 
to a = — 1 (i.e. riT = 0). In general the power law spectrum just assumes the absence of features in the 
spectrum of gravitational waves at the wavelengths of our interest. 

In practice, when considering pulsar timing, we are interested in calculating the expected mean square 
deviation of the timing residuals due to stochastic background of gravitational waves (jl9bp . In order to 
evaluate < R^{t) > from expression (jl9bp we require to specify the limits of integration k^in and k^ax- 
kmin and kjnax determine the frequency range of gravitational waves that can be probed by pulsar timing 
measurements. The lower limit kmin is determined by the time duration of observations Tobs, kmin ~ 
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2'Kfobs/c = 2'K/cTobs- In our estimates we shall assume Tobs ~ 10 yrs. The upper limit kmax ~ ^ixjcbt is 
determined by the duration of single observation it (in other words, the time of integration), which is usually 
of the order of 1-2 hours. We note here that it is this time (and not the time between consecutive observations, 
of order of weeks) which determines k^ax in timing residuals. Indeed, if the period of a gravitational wave 
is smaller than (5t, its effect is smeared out by averaging procedure. But if the period of the gravitational 
wave lies between the averaging time and the sampling time, the wave will clearly manifest itscft in the 
timing residuals. Some authors erroneously use the inverse sampling time as kmax, apparently guided by 
the analogy with time series analysis. Thus, in our case, it is safe to assume 6t <^ Tobs (i-c kmax 3> kmin), 
and set kmax ~ oo in numerical evaluations below. Furthermore, we shall be working under the assumption 
kD ~ 2-KD/Xgw S> 1, which corresponds to the reasonable assumption that the gravitational waves of our 
interest {\gw ^ 10 ^yrs) have wavelengths much shorter than the distance to the pulsar {D ^ 10 kpc). 

As can be seen from expression (|23p and the considerations in Appendix the behaviour of the transfer 
function R^{k) depends on value of the quantity 'SiTe^kD/2. In order to analyze the various possibilities let 
us introduce 

-1/2 



3.2x10 



-3 



10 kpc 
D 



Tob 



10 yrs 



(28) 



Below we shall analyze the two possibilities, e ^ and e 3> e*, separately. 

In the case e <C in transfer function i?^(fc), in expression (j23p . we can neglect the second term in the 
square brackets in comparison with the first. Furthermore, in the term sin^ [kcT [1 — e)) /2 we can neglect 
the rapid oscillatory factor. Thus, for the transfer function we get 



i?2(fc) 



3fc2c2 



1 - cos(fccr(l - e)) 



-T^e^kD 
2 



3k^c^ 



(29) 



Substituting the above approximation (j29p . taking into account the definition (|24p and a power law spectrum 
(|26p . into expression (|19bp . and setting the limits of integration as mentioned above, we arrive at 

Tobsf^lifobs) 



< R\t) >f 



for e ^ e*. 



(30) 



247r2 (1 - a) 

In the case e ^ e*, neglecting the first term in the square brackets with respect to the second in 
and ignoring the rapid oscillatory factor, the transfer function can be approximated as 



i?2(fc) 



3fc2c2 



1 - cos(A:cT(l - e)) 



1 



Tie^kD 



2fcc2 ■ 



In this case, the expression for (|19bp takes the form 



< R\t) 



127r2(l-2a) \e 



for e » 



(31) 



(32) 



Comparing expressions ([30)) and ([32]) it can be seen that when e ':$> e^, the surfing effect leads to a strong 
resonance contribution (proportional to kD) in the timing residual compared with the case when e <C e*. 
This dominant resonance contribution comes from gravitational waves traveling at an angle cosO ~ (1 — e) 
to the direction of signal propagation from the pulsar (sec Appendix |X] for details). 
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From expressions pop and (|32p . it follows that the direct measurement of pulsar timing residuals would 
be able to measure or constrain either he or hce, depending on the value of e compared with e^. A null result 
in timing residual measurements would place the following upper limits 



he < 4.9 X lO"^-' 



Rr 



\ flO yrs 



0.1 /iscc 



To. 



bs 



for e ^ e* 



or 



hrC < 1.1 X 10" 



Vl -2a 



Rr 



0.1 ^sec 



^ f 10 yrs 



10 kpc 
D 



for e » e. 



(33) 



(34) 



where Rrms = \/<RHt)> is the precision of the pulsar residual timing, and he = hc{fobs) is evaluated at 
fobs = 0.1 yrs^^. It is also convenient to present this limits in terms of the density parameter fig 



''gw 



^gw < 5.3x10" 



(1 - nr/A) 



Rr 



f 10 yrs 



or 



fig^e^ < 4.0x10" 



(1 - «t/3) 



0.1 /ISCC / \ Tobs 

10 kpc 



for e <C e* 



D 



Rrms 

0.1 /isec 



^ 10 yrs^ ^ 



To, 



■6s 



for e e* 



(35) 



(36) 



Thus, from (|33|) (or (j35p '). it can be seen that for e <C e*, when surfing effect is not important, pulsar timing 
sets limits directly on he (or equivalently on f2gu)), i-e. the strength of the gravitational wave background. On 
the other hand, when €':$> t^, and surfing effect becomes dominant, from (or ([55]) ). it follows that pulsar 



timing sets limits on the combination hct (or ^.g^e^ equivalently). The upper limits from pulsar timing, 
along with possible sources and sensitivity levels of various experimental techniques to detect gravitational 
waves, are illustrated on figure [H 

As follows from the above discussion, and can be seen from figure [21 an independent knowledge of he 
would enable us to directly constrain the parameter e, i.e. constrain the deviation of speed of gravitational 
waves from speed of light. From expression (|34p we arrive at the following constrain on e 



e < 1.1x10" 



Vl - 2a 



/lO 



-15 



10 kpc 
D 



Rrms \ / 10 yrs 



0.1 ^sec 



To, 



■bs 



In terms of the density parameter f2g^ the constraint has the form 



e < 6.4x10" 



10 



-10 



gw 



/lOkpcy / Rrms 

\ D ) VO.l Afscc 



10 yrs 

To. 



bs 



(37) 



(38) 



In the next section we shall discuss the various viable candidates for a stochastic gravitational wave 
background and explicitly calculate the achievable limits on e. We shall also discuss the implications of the 
surfing effect for theories with massive gravitons. 



V. THE PHYSICAL IMPLICATIONS OF THE SURFING EFFECT 

The analysis in Section IIVI indicates that the surfing effect in pulsar timing can yield interesting con- 
straints on e parameter and consequently the mass of graviton in a sufficiently strong gravitational wave 
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FIG. 2: The upper limit on strain amplitude he and velocity parameter e for gravitational waves, achievable 
by pulsar timing residual measurements with precision Rmis = 0.1 fiscc and time of observation Tots = 
10 yrs. The shaded area shows the region that can be probed or ruled out by pulsar timing observations. 
The horizontal lines show the strain he, at / = (10 yrs)^^, for some viable sources of gravitational wave. 

background with flguj ^ 10^^" (see ([55]) '). It is important to note that this method is fundamentally limited 
by the value e*, which is currently about 3 x 10^"^ (see ^28\i ). Although an increase in the time of observation 
will improve overall precision, it will also increase the value of e^, thus worsening the potential constraints 
on e. In future, the method can become more sensitive with implementation of large radio telescopes like the 
Square Kilometer Array (SKA) (see [l^ for detailed discussion of SKA and its usage in pulsar astrophysics), 
which would improve the limitations to e* ~ 10^'^. Furthermore, as seen from expression (|38p (or (|37p ). 



increasing the pulsar timing accuracy (for example, using pulsar timing ensembles 4l|) can reduce the limit 
down to the critical value e*. 

The gravitational wave background, at the frequency range of our interest {fgw < O.lyrs"^), consists of 
contribution from a variety of well established astrophysical and cosmological sources Q as well as possible 
contribution from exotic remnants of early universe [3| , [3| ■ The strongest contribution to the gravitational 
wave background, at these frequencies, come from the background of extragalactic coalescing supermassive 
binary black holes (SMBH) Q, Q, H- For this reason, below in the subsection IV Al we shall study 
the implications of the surfing effect for this background. Following this, in subsection I VB|, we shall analyze 
the consequences of the limitations on e for theories with massive gravitons. 
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A. Gravitational wave background from extragalactic black holes 



As was mentioned above, one of the strongest sources for a stochastic gravitational wave background at 
frequency range of our interest, fg^ ^ T^obs ~ yi's^^, conies from the extragalactic black hole binaries. 
Various groups have conducted a theoretical study on the strength of this background [l^, [l^], 
There is a general consensus on the expected gravitational wave strain for this background 

_ 2 

hc{f)^lQ-''{-^\ \ (39) 



l^Hz 



corresponding to the value for the density parameter 



17g^(/)«2.4xl0-^"(— . (40) 



^0.1 yrs-i 

The uncertainty surrounding this value of he arises mainly due to the uncertainty in the galaxy merger rates 
as well as some other astrophysical factors. Taking into account these uncertainties, the amplitude lies in 
the interval h^f = l^iHz) w 2.5 x 10"" - 4 x 10"!^ Q- 

The expected strain he from the background of SMBH allows to place significant bounds on the e 
parameter. Substituting expression (|39p into expression (j37p . and setting a = —2/3, we arrive at the 
following limit on upper e 

10 kpc\ 2 / Rrins \ / 10 yrs \ 2 



e < 3.7x10^3 



(41) 



D J \0.1 /isccy \ Tabs 

Thus, the stochastic gravitational wave background of extragalactic SMBH mergers can potentially place 
very stringent constraints, e < 0.4%, on the speed of gravitational waves. 



B. Implications for theories with massive gravitons 



The phenomenological parameter e is directly related to the mass of the graviton nig (see ([6])). It is 
convenient to rewrite expression ^ in the form 



e(fc) = eo 



where to 



k J' 2nh 
For the fiducial strength of gravitational wave background we get 

10-1" \^ /ID kpc\^ / R 



eo < 8.3x10 



-3 



Vl - nr/B 



D 



W 10 yrs 
0.1 /iscc / 



(42) 



(43) 



Note that the factor tit/^, in the above expression (|43p (compared with the factor rix/i in expression ([55]) ). 
arises because we are constraining (compared with constraints on e in (j38p ). This leads to an extra factor 
(k/kmin) in integral (jl9bp and hence slightly modifies the result. The above limit on Eq implies the following 
limit on the mass of the graviton 

'10-i"\ VlO kpcV ^ Rrms WlOyrsV 



rug < 1.1x10 



-25 



eV 



^1 - nT/5 



gw 



D 



0.1 Msec ) \ To, 



obs 



(44) 
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From expression (j42p it follows that, stronger constraint on trig require smaller values of kmin, i-e. require 
a longer time of observation Tots- On the other hand, the strongest possible constraint for is determined 
by the value of (which increases with the time of observation, see expression (j28p '). For this reason, an 
increase in Tots beyond a value of approximately 25 yrs will not lead to an improvement in constraining nig. 

As a concrete example, let us assume that the gravitational wave background from SMBH coalesces 
dominates at frequencies 0.1 — 1 yrs^^, and that its properties are not affected by the non-zero mass of 
graviton. Then the existing four-years precise timing of PSR B1937-I-21 4l| allow to significantly constrain 



the mass of the graviton. Setting Tohs ~ 4 yrs, Rrms — 0.17 fiscc, D = 8.3 kpc, ut ~ 2/3 and f7g.iu(T^ 
4.2 X 10^^" (see pO)) ) in expression (|44)). we arrive at a limit 

rUg < 3.6 X 10^2^ eV, 



bs' 



(45) 



corresponding to a Compton length for graviton of Ag 



h/mgC > 3.4 X 10"'^^ km. This bound is three orders 



42l | and is comparable to future limits from SMBH 



stronger than the current limit from Solar system tests 
mergers obtainable with LISA (see 43| and references therein). It is worth stressing, that the limits from 
pulsar timing are more robust and less model dependent than the prospects for LISA. 

The surfing effect in pulsar timing puts stringent constraint on the mass of graviton in some theories of 
gravity (see [4j| ) . In 45| the authors propose massive gravitons as a viable candidates for cold dark matter 
in the galactic halo. At the frequency ranges of our interest, these massive gravitons imply e « 0.5. The 
existing precise timing of PSR B1937+21 place direct limits on the parameter flgw^^ ^ 2x10"^"^ (setting 
Rrms ~ 0.17 /iscc and Tots ~ 4 yrs in expression (|36p ) . This implies that massive gravitons, as candidates 
to explain the dark matter in the galactic halo, can be ruled out with the current observations. 



VI. CONCLUSIONS 



In this work we have analyzed the consequences of the surfing effect, introduced in |30[, for pulsar timing 
observations. The surfing effect, due to the transverse nature of gravitational waves, leads to a strong 
observable signature only when the speed of gravitational waves is smaller than the speed of light. In order 
to analyze this possibility, we have introduced a parameter e, which characterizes the deviation of speed of 
gravitational waves from speed of light. By studying the pulsar timing residuals in the presence of a single 
plane monochromatic gravitational wave, followed by a generalization to an arbitrary gravitational wave 
field, we show the presence and importance of surfing effect in the case when e =/= 0. 

The surfing effect allows to place significant bounds on the parameter e. For a timing accuracy of 
Rrms ~ 0.1 yusec, and assuming a realistic background of gravitational waves from extragalactic super 
massive black hole binary mergers, the achievable limits are e < 0.4%. The strongest achievable bounds on 
e are determined by e*. For a pulsar at a typical distance D = 10 kpc the value is w 0.3%. This limit 
could potentially be slightly improved by observing pulsars at a greater distance D. 

The surfing effect leads to interesting consequences for theories with massive gravitons. Using the 
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existing observations, we have constrained the mass of graviton to rUg ^ 4 x 10~^^ eV, which is three 
orders of magnitude stronger than the current limits from Solar system tests. With future observations this 
constraint could improve by an order of magnitude. Based on the existing observations, we have also ruled 
out massive gravitons as candidates to explain the dark matter in the galactic halo. 

In comparison with precision intcrferomctry methods considered in jsoj l. pulsar timing measurements 
(due to their high precision) should be able to put tighter constraints on e. In any case, these two methods 
of constraining e are independent and hence should be considered complementary. 
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APPENDIX A: EVALUATION OF THE TRANSFER FUNCTION 



Let us evaluate the integral in expression ^2. 

+1 



I(k)= / 



sm 



(1 - e - 



(Al) 



in the physically interesting case when e ^ and kD ^ 1. The integral can be separated into two distinctive 
contributions 



I{k) = lNR{k)+lR{k), 



(A2) 



where iNnik) is the non-resonance contribution 



(fc) = 



(l - 



+1 



dfi (l — fi 



(A3) 



and Aa|j(fc) is the resonance (or, in other words "surfing",) contribution 



(A4) 



Aa|(fc) = I df, (l-/i2)^ 

The quantity Afi occurring in the limits of integration in the above expressions is fixed by the condition for 
the resonance to occur. This condition corresponds to the region, around /i = 1 — e, where the sine function 
undergoes a few oscillations. Thus Afj, = NXgw/D = 2'KN/kD, where N is the number of oscillations of the 
sine function, around the point = 1 — e, included in evaluation of the resonance. The value of N is limited 
by the condition = 2nN/kD <^ e, implying N <C ekD /2Tr. Since in all our considerations we assume 
e ^ 1, and e^kD ^ 1, the condition imposed on N is consistent with an additional condition ^ 1 that 
we shall assume. 

When evaluating (jASp . since we assume e ^ 1, we can neglect the second integral in comparison with 
the first. In evaluation fo the remaining integral we can set e = 0. Thus, we get 

1 



(k) 



-1 

1 

^ fdfi (l + /i)' (l- cos(fci?(l-M))) 



dn (1 + /i) 



(A5) 



where, assuming kD 3> 1, we have explicitly separate out the rapid oscillatory part and neglected it in the 
last line. 
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In order to evaluate (|A4[) . in the case of e ^ and kD 3> 1, it is helpful to notice that the factor 
(l — /i^)^ in the right side of (|A4p is a slowly varying function over the range of integration. Taking this 
factor (evaluated at /i = 1 — e) outside the integral we get the following approximation for the resonance 
part of the transfer function 



dfi 



sm 



2TTe^kD {l-O 



1 



« 27re^kD 



= 2e^kD 



dx 



■ 2 

sm X 



-N-K 



(A6) 



Finally, the total transfer function, given by the sum of the non-resonance (|A5P and resonance parts 
(|A6p . has the following form 



I{k) = lNR{k)+lR{k) 



3 , 



(A7) 
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